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Abstract 

The nonlinear dynamics of a bundle of magnetic flux ropes driven by sta- 
tionary fluid motions at their endpoints is studied, by performing numerical 
simulations of the magnetohydrodynamic (MHD) equations. The develop- 
ment of MHD turbulence is shown, where the system reaches a state that is 
characterized by the ratio between the Alfven time (the time for incompress- 
ible MHD waves to travel along the field lines) and the convective time scale 
of the driving motions. This ratio of time scales determines the energy spectra 
and the relaxation toward different regimes ranging from weak to strong tur- 
bulence. A connection is made with phenomenological theories for the energy 
spectra in MHD turbulence. 
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I. INTRODUCTION 



The search for universal scaling properties in magnetohydrodynamic turbulence, of the 
type developed by Kolmogorov [1] for hydrodynamic turbulence, has been the subject 
of many theoretical and numerical studies, since the pioneering phenomenological argu- 
ments put forward independently by Iroshnikov and Kraichnan [2,3] (hereafter IK). The 
phenomenological arguments of Kolomogorov, which give rise to the famous power law 
Ek ~ k~ 5 ^ 3 for the omnidirectional energy spectrum in terms of the wavenumber k, were 
modified by IK to include magnetic field effects, deriving an energy spectrum power law 
like Ek ~ A;~ 3 / 2 . An alternative point of view, however, was proposed [4] indicating that 
the original Kolmogorov scenario is still applicable in certain cases to MHD as well as to 
hydrodynamics. Along similar lines, a phenomenological theory was proposed in Ref. [5] for 
a steady inertial range spectral law that reduces to the IK and Kolmogorov laws in appropri- 
ate limits. To distinguish between the two power laws is a difficult task both observationally 
and numerically [6], partly because the power law indexes are too close, but also because 
assumptions on the theories, such as homogeneity, isotropy and time stationarity are often 
lost. 

For instance, when a strong, externally supported, magnetic field is present, the key 
assumption of isotropy breaks down. It has been shown [7-10] that the cascade renders itself 
anisotropic, and that the spectrum in the direction of the magnetic field becomes strongly 
suppressed. This situation also holds for compressible MHD, where a variance anisotropy 
is also found [11], i.e., components of the fluctuating fields in the strong magnetic field 
direction are small compared to the transverse components. On theoretical grounds, this 
configuration has been found [12] to be appropriately described, in the low-frequency limit, 
by the so called Reduced MHD [13,14] approximation (RMHD). 

Phenomenological theories alternative to the IK scenario [2,3] when a strong magnetic 
field is present have been also presented to describe the perpendicular spectrum in anisotropic 
conditions [15,9]. Recently, a formal closure model has been introduced, which obtains a 
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kinetic equation for the anisotropic energy spectrum in the limit known as weak MHD 
turbulence [16]. 

Besides the anisotropy introduced by a strong magnetic field, another effect that may 
modify the cascade situation is the presence of driving boundaries. We intend to address 
the effect of both a strong external magnetic field and the driving boundaries on MHD 
turbulence, by considering the particular case of a plasma under a strong uniform field in 
the z-direction, limited within the transverse planes z = and z = L z where an imposed 
velocity field is applied. This case is inspired in the theoretical model of a coronal magnetic 
loop, which is forced through convective motions at its footpoints. Such a model loop, 
proposed in Ref. [17], has been widely studied in many different contexts, although not 
necessarily in a turbulent scenario. Most studies address the issue of current sheet formation 
[18-20], non-steady reconnection [21,22] and coronal heating [23,24,21,25,26]. However, this 
paper is not aimed at studying the dynamics of coronal magnetic loops. Instead, we focus 
our analysis on the (steady) perpendicular energy spectrum of boundary driven RMHD 
and make a connection with the phenomenological theories of MHD turbulence. We pay 
particular attention to the issues of numerical resolution required to yield well resolved 
turbulent spectra from direct numerical simulations. 

The organization of the paper is as follows: in section II we present our model and the 
RMHD equations. Section III contains a description of the boundary conditions and the 
numerical code. Section IV describes the statistically steady regime and the dissipation rate 
scaling law. In Section V we show the numerically obtained energy spectra and discuss the 
connection with phenomenological theories of turbulence. Issues of numerical resolution are 
also addressed. Section VI discusses the typical dissipative structures. Section VII contains 
the conclusions. 
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II. MODEL AND EQUATIONS 



Let us assume a low-/? magnetofluid (/? is the ratio of gas pressure to magnetic pressure) 
permeated by an initially uniform magnetic field B = B z, which is elongated along the 
z-direction (i.e. L z ^> L±), as shown in Figure 1. Under these conditions, the so-called 
reduced magnetohydrodynamic equations (RMHD) are applicable [13,14,12] to describe the 
dominant low-frequency non-linear motions. Within this approximation, the magnetic and 
velocity fields are both divergence-free and can be expressed as 

B = B z + V ± x (az) (1) 
u = V x x (ipz) (2) 

where a(x, y, z, t) is the magnetic flux function and ip(x, y, z, t) is the stream function. The 
RMHD equations in terms of these scalar potentials are: 

d t a = v A d z ip + [ip, a] + rj'V ± a (3) 

d t w = v A d z j + [ip, w] - [a, j] + vV\w (4) 

where w(x,y, z,t) = —V^ip is the parallel vorticity and j(x,y,z,t) = — V\a is the parallel 
current density. The coefficient v A = Bq/^Attp is the Alfven speed, v is the kinematic 
viscosity and r\ is the magnetic diffusivity. The non-linear terms in these equations are 
expressed in terms of standard Poisson brackets, i.e. [u,v] = (d x u)(d y v) — (d y u)(d x v). 

This particular theoretical setup is relevant to several plasma applications, ranging from 
tokamaks to magnetic loops in the solar corona. The RMHD framework has been used to 
study the dynamics of coronal loops in solar active regions [23,19,21,25]. The axial (and 
approximately constant) magnetic field in these loops has both ends or footpoints anchored 
in the solar photosphere. The photosphere is a high-/? plasma, which is also convectively 
turbulent. These convective motions at the photosphere, in turn move the magnetic fieldlines 
around, and drive the coronal part of the loops into a rather complex dynamical scenario. 
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We specify the velocity fields at the boundaries as 



i>{z = 0) = 0, 



i>{z = L z ) = *(x,y) 



(5) 



where ^(x, y) is the stream function which describes stationary and incompressible footpoint 
motions. The strength of this external velocity field is proportional to a typical velocity U p . 

To transform equations (3)- (4) into their dimensionless form, we choose l p = L±/(2ii) 
and L z as the units for transverse and longitudinal distances. Since the dimensions of all 
physical quantities involved in these equations can be expressed as combinations of length 
and time, let us choose t A = L z /v A as the time unit. The dimensionless RMHD equations 
are: 



where S 1 = %r and R 1 = ^ are the (dimensionless) magnetic and kinetic dissipation 
coefficients. Hereafter, we will consider the case S = R, and thus the (common) dissipation 
coefficient will be the only dimensionless parameter explicitly present in equations (6)- (7) . 

III. DESCRIPTION OF THE CODE AND BOUNDARY DRIVING 

We numerically integrated equations (6)-(7). To this end, ip and a are expanded in 
Fourier modes in each (x,y) plane (0 < x, y < 2n and < z < 1). The corresponding 
Fourier coefficients ipk(z, t) and a^(z,t) are evolved in time using a semi-implicit scheme: 
linear terms are treated in a fully implicit fashion, while nonlinear terms are evolved using 
a second order Runge-Kutta scheme. Also, nonlinear terms are evaluated following a 2/3 
fully dealiased (see Ref. [27]) pseudo-spectral technique. 

To compute ^-derivatives we use a standard method of finite differences in a staggered 
regular grid (see for instance Refs. [13,19]) of N z + 1 points. The stream function is computed 



d t a = d z ip + [ip, a] + - V 2 ± a 



(6) 



d t w = d z j + [ip, w] - [a,j] + —V\w 



(7) 
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on points Zi = i/N z (i = 0,...,N Z ), while the magnetic flux function is computed on 
z i+i/2 — (i + 1/2) /A^ (i = 0, . . . , N z — 1). Boundary conditions for the stream function ip 
are given at the plates z = and z — 1. Therefore, equation (7) is not integrated on these 
planes, but it is evolved in time in all the internal gridpoints Zi = i/N z (i — 1, . . . , N z — 1). 
We specify the stream function in equation (5) as 

*k = *o = W*i» if 3<kl p <4 (8) 
^ k = elsewhere 

This expression imitates a stationary pattern of eddy motions of diameters between 
L_[_/4 and Lj_/3, rotation speeds U p and typical turnover times t p = l p /U p . Our choice of 
this narrowband and non-random forcing ensures that the broadband energy spectra that 
we obtain are exclusively determined by the nonlinear nature of the MHD equations. 

The strength of the external driver is quantified by the dimensionless factor \1> = t A /t p in 
equation (9), which is given by the ratio of the Alfven time of the system (i.e. the response 
time to an impulse applied at the boundary) to the timescale of the driver itself (i.e. the 
eddy turnover time). Note that the velocity fields that we are imposing at the boundaries 
are stationary. For non-stationary footpoint motions, which might for instance represent 
wave activity or time-correlated random flows, other timescales should be considered in 
connection with the driver. 

IV. STATIONARY REGIMES AND SCALING LAW 

As mentioned in Section 2 , the dissipation coefficient S (with S — R) is the only 
dimensionless parameter present in equations (6)- (7). Just as important is the dimensionless 
factor t A /t p introduced by the external force applied at the boundary. Therefore, we are left 
with these two dimensionless numbers to characterize the solutions of the RMHD equations 
[i.e. equations (6)-(7) with the boundary condition given by equation (5) and equation (9)]. 

From purely dimensional considerations, we know that for any physical quantity, its 
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dimensionless version Q should be an arbitrary function of the only two dimensionless pa- 
rameters of the problem, i.e. 

Q = ^(Qi,Q 2 ), Qi = j-, Q2 = S (9) 

For instance, for the dissipation rate per unit mass e 

Q = 4 = ^,S) (10) 

A sequence of numerical simulations for different values of the dissipation coefficient 
performed in Ref. [25] (see also Ref. [28]), indicated that the dependence of the dissipation 
rate with S is rather mild. This relative insensitivity of the dissipation rate with the dis- 
sipation coefficient is consistent with similar results obtained from experiments in purely 
hydrodynamic turbulence [29-31]. It is also consistent with one of Kolmogorov's hypothesis 
for statistically steady turbulent regimes at very large Reynolds numbers, which assumes 
that the dissipation rate remains finite in the limit of vanishing viscosity [31]. Therefore, 
if we assume that the dissipation rate in statistically steady turbulent regimes in RMHD 
is approximately independent of the dissipation coefficient S, we readily obtain that the 
dissipation rate will only depend on the time ratio t A /t p . A second series of simulations 
[25], in this case for different values of the ratio t A /t p , led to the following expression for the 
function J 7 , 

e=|(^) S , s = 1.51 ±0.04 (11) 

In the present paper, we performed numerical simulations of equations (6)- (7) with S = 
2000 and different values of the ratio t A /t p in the range [0.1, 1]. Note that t A /t p can also 
be written as t A /t p = (L z /l p )/(v A /U p ). For consistency with the RMHD approximation, 
both L z /l p and v A /U p should remain much larger than unity. However, as mentioned above, 
the only relevant parameter in this problem is t A /t p , which is free to take any value. An 
additional simulation with t A /t p = 2 and S = 800 is also performed. Numerical resolution 
ranges from 256 x 256 x 64 to 2048 x 2048 x 128 gridpoints. A Beowulf PC cluster is employed 



to perform those runs, and a parallel RMHD code has been designed. Since finite-differences 
are used in the z-direction, while a pseudospectral method is employed for the transverse 
directions, a very efficient parallelization is achieved by performing the transverse gradients 
locally in each machine, while employing communication between machines only to perform 
the ^-derivatives. The typical behavior of the magnetic and kinetic energy as a function of 
time is shown in the top panel of Figure 2, while the bottom panel shows the dissipation 
rate for the particular case of t A /t p = 0.5 with 256 x 256 x 64 resolution. After an initial 
transient, the dissipation rate is seen to approach a statistically steady level. 

V. ENERGY SPECTRA 

Energy spectra of turbulent fluctuations are essential both for phenomenological and 
statistical theories of MHD turbulence. We compute the perpendicular spectra at each 
plane and integrate in z to obtain a perpendicular spectrum, 

Ek ± = lf E [\h(k x ,k y ,z)\ 2 + \u(k x ,k y ,z)\ 2 ]dz (12) 

where the total energy can be obtained as E = J E k± dk±. 

As we have pointed out in the previous section, the system reaches an approximately 
steady state where the forcing is compensated by dissipation (Figure 2 bottom). The first 
stage of the system evolution is followed using a low resolution simulation. Once the steady 
state is achieved, a higher resolution simulation is started, by taking the Fourier coefficients 
from the previous run and padding the additional high k coefficients to zero. After a short 
transient, those additional coefficients are populated by the nonlinear cascade in a few 
eddy turnover times. In this way, a high resolution steady state spectrum is obtained in a 
relatively shorter computer time. The high resolution runs (2048 x 2048 x 64) are however 
computationally demanding, even for runs of a few eddy turnover times. 

In Figure 3 we show the magnetic and kinetic energy perpendicular spectra for different 
values of the quantity t A /t p . The right panel plots show the compensated spectra (magnetic 
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plus kinetic), which allow determination of the power law index a observed in the inertial 
range where E k± ~ k^_ a . A first clear result observed in these plots is that the power law 
index depends on the value of tA/t p . For low values of t A /t p a rather steep spectra (a > 2) 
is observed. At intermediate values of tA/t p the power law index is approximately equal to 
a = 2. For larger values of t A /t p the spectral law becomes flatter (a < 2). As shown in 
Figure 4, when t A /t p = 2 a value close to a = 5/3 is obtained. For this last simulation a 
lower value of S — 800 has to be employed, in order to have a well resolved spectrum (see 
discussion below). 

Phenomenological theories [15,32] and a more formal theory based on closures [16] have 
been proposed to explain a k~ 2 spectra in a regime known as the weak turbulence limit. An 
alternative phenomenological derivation of this power law can be done using the general 
approach proposed in Ref. [5] in terms of the triple correlation time. The energy supply e 
is assumed to be equal to the /c-independent spectral transfer rate. As stated in Ref. [5], if 
r 3 (k) is the time scale for decay of triple correlations, which induce the spectral transfer to 
wavenumbers higher than k, then [3,5] 

kE k 

*~«l>) m (13) 

where r n i(k) ~ (kuk)' 1 is the non-linear time. Equation (13) implies a spectral transfer 
time t s ~ t 2 Jt 3 . If r 3 (k) ~ T n i(k) and equipartition between magnetic and kinetic energy 
is assumed E k ~ E k ~ u 2 k /k, the usual E k ~ k~ 5 ^ 3 Kolmogorov scaling is recovered. When 
r 3 (k) ~ T A {k) = (/c6) _1 , where b is the rms of magnetic fluctuations, and isotropy and 
equipartition is assumed, the IK scaling E k ~ k~ 3 ^ 2 is obtained. If there is a strong external 
magnetic field the energy cascade is anisotropic and r nl ~ (kx_u kA )~ l . For this case, the 
appropriate triple correlation time r 3 is the Alfven wave crossing time (in the direction of 
the external field) r 3 ~ t A = L Z /V A , a fixed external quantity independent of k± (this is 
different from the /c-dependent Alfven time r A {k) associated with fluctuations in an isotropic 
situation). Again assuming equipartition, 

e ~ t A k ± E k± k 2 ± u 2 k± ~ t A k 4 ± E 2 k± (14) 
9 



and a E k± ~ kj 2 spectra is obtained. This is equivalent to the wave packet interaction 
approach presented in Ref. [15], based on the IK scenario [2,3] but under strong anisotropic 
conditions. It is also interesting that an essentially similar argument [33] has been applied 
to the case of purely hydrodynamic turbulence in the presence of strong rotation, where 
the triple correlation time is set equal to a fixed external timescale T3 ~ tn = 1/fl This 
suggests an analogy between strong rotating turbulence and RMHD (i.e. MHD turbulence 
with a strong mean field). The simulations presented here show that the weak limit regime 
of E k± ~ k~ 2 is obtained when ^/i p « 1/2. For larger values of this timescale ratio, the 
energy spectrum approaches a Kolmogorov-like regime E k± ~ k ± 5 ^ 3 . 

For the small tA/t p = 0.1 case the spectrum is E k± ~ k~ 2A which would indicate that 
the weak turbulence regime gave rise to a regime in which the perpendicular cascade is 
suppressed even more. In these numerical experiments we have kept the value of S = l 2 /ij]t A ) 
fixed, while changing the value of t A /t p = (L Z U P ) /{l p v A ). If we consider that all parameters 
are fixed, except for U p , we can interpret these runs as different evolutions of the system 
when varying the intensity U p of the boundary driving. The result is that as the driving 
intensity is enhanced (i.e. larger tA/t p ), the spectrum becomes flatter and the regime goes 
through weak turbulence to a strong Kolmogorov-like state. As shown by Figure 3 (left 
panels), the kinetic energy spectrum is always smaller than the magnetic energy spectrum 
and the ratio E k /E k is smaller (especially at the low k structures) when the forcing is weaker. 

The phenomenological argument expressed by equation (13) can be modified to address 
the very weak forcing case. For that case, the energy at the low and intermediate k values is 
essentially magnetic E k± E k± (see Figure 3). For the non-linear time r n i(k) ~ (kuk)' 1 , we 
replace the /c-dependent velocity u k by a constant which we expect to be of the order of the 
imposed boundary weak velocity U p (i.e. r n i ~ ik^Up)^ 1 ). A scale invariant characteristic 
speed corresponds to a kinetic energy spectrum E k± ~ kj 1 . Again assuming r 3 ~ in 
equation (13) we find, 

e ~ t A k ± E b k± k 2 ± U 2 (15) 
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and a steep spectrum E% ± ~ A; J 3 is obtained. We should note that there are theoretical 
developments predicting spectra steeper than k~ 2 , due to nonlocal effects [15], or due to 
high values of the correlation between magnetic and velocity fields [34]. The latter one 
corresponds to (u 2 ) ~ (6 2 ), a situation that would not be applicable for the present case. 
The effects of the finite values of S considered here are also expected to be more important 
for very weak forcing. 

Phenomenologies [15] and closures [16] are also usually presented in terms of the Elsasser 
fields z ± = u±b, so for completeness we computed the corresponding perpendicular spectra 

E t = -j E \^{k X ,ky,Z)\ 2 dZ (16) 

fc ± <fc2+fcg<fc ± +Afc 

The results for the cases tA/t p = 0.1, 0.5, 1 are shown in Figure 5. The compensated spectra 
for the product are shown on the right panels of Figure 5. Weak turbulence theory 

[16] predicts an index a — 4, which for this particular problem corresponds again approxi- 
mately to the case t^Ap = 1/2- A higher power index a ~ 4.9 is instead obtained for the 
very weak forcing case t A /t p = 0.1. 

A rather technical but important issue when determining spectra is whether the dissi- 
pation range is well resolved. We have considered here an exact laplacian dissipation term 
in the dynamical equations (i.e. no hyperviscosity is used) which makes harder to achieve 
higher spatial resolutions. An example from a simulation with t^/tp = 1 is shown in Fig- 
ure 6. A /c-dependent Reynolds number is computed as R k = u k /{k u) = ^ (E^ / k) / v . The 
dissipative region in /c-space is characterized by R k ~ 1. The fact that R k reduces mono- 
tonically and falls below unity at the largest values of k, is compelling evidence that the 
dissipation range is well resolved in our simulations. Similar results are obtained for the 
different values of i^/i p , although higher perpendicular resolution is required as this value is 
increased. For t^/tp = 2 the excessive computational requirement has forced us to consider 
a lower value of S — 800 to resolve the dissipation range. Another direct resolution study is 
shown in Figure 7 where different perpendicular spatial resolutions have been adopted for a 
simulation with t^Ap = 0.1. What this figure shows is that the spectrum is better resolved 
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and extends to higher wavenumbers as the resolution is increased. However, a very good 
picture of the lower k part of the spectra can still be obtained with the low resolution runs. 

VI. SPATIAL STRUCTURES 

Figure 8 is a cross section of the simulation box (perpendicular to the strong external 
magnetic field) showing the spatial distribution of electric current density (z-component) at 
t = 120t A for the cases t A /t p = 0.1 (Figure 8a) and t A /t p = 0.5 (Figure 8b). Intense positive 
z currents correspond to white regions, while intense negative z current concentrations are 
indicated in black. The current density structures extend in the z-direction (not shown) 
with almost no variation (as consistent with the RMHD approximation) in a current sheet- 
like form. These structures are highly dynamic in the x, y plane and evolve in time. They 
indicate the presence of multiple reconnection events (where the reconnecting magnetic field 
lies in the plane perpendicular to the external field) occurring within the turbulent non-linear 
dynamics. 

The minimum thickness that current sheets can reach is determined by the magnetic 
diffusivity. In numerical simulations of turbulent regimes, the magnetic diffusivity is usually 
made as small as possible, in such a way that the sheet thickness (which is about the smallest 
spatial feature expected in these simulations) is just marginally resolved. The two cases 
presented in Figure 8 share the same value of magnetic diffusivity (S = 2000), and in fact 
it is possible to find very narrow current sheets in both cases, close to the resolution limit. 
Both cases correspond to statistically steady regimes. Nonetheless, the case corresponding 
to the less intense forcing (t A /t p = 0.1 in Figure 8a), shows a comparatively fewer number 
of current sheets. Therefore, the number of current sheets seems to be controlled by the 
intensity of the forcing. The widths of these current sheets on the other hand, are mostly 
determined by the dynamics. As a result, in both cases a wide range of current sheet widths 
can be observed. 
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VII. CONCLUSIONS 



We presented results of the turbulent dynamics of a bundle of magnetic flux ropes driven 
at their endpoints with steady convective motions. Two dimensionless parameters control- 
ling the response of the system are identified: the dissipation coefficient and the ratio between 
the Alfven wave box crossing time to the forcing timescale, t A /t p . As we showed in a previous 
work [25], the dissipation rate strongly depends on the value of this timescale ratio, and is 
only weakly dependent on the dissipation coefficient. It is found that the external driving 
generates a broadband perpendicular energy spectrum with a slope which is determined by 
the timescale ratio t A /t p . For small values of this ratio, a spectrum with slope a —2.4 is 
obtained. For t A /t p ~ 0.5 the spectrum is ~ kj 2 , which would indicate a realization of the 
regime known as weak turbulence. For higher values of t A /t p the spectrum approaches the 
Kolmogorov form k^ 5 ^ 3 . Both the weak turbulent and the Kolmogorov spectra can be ob- 
tained from a common phenomenological framework based on the assumed value of the time 
for decay of triple correlations. The regime corresponding to spectral power laws with slope 
-2.4 and lower for low values of t A /t p should be attributed to the relatively small amplitude 
flow at the boundaries. The phenomenological framework can be modified to consider this 
weak velocity field case and a steep magnetic spectrum ~ kj 3 is obtained. We also show the 
development of small scale spatial structures, in the form of current sheets oriented along 
the mean field. The number of current sheets seems to be controlled by the timescale ratio 
tA/tp, while the widths of these current sheets on the other hand, are mostly determined by 
the dynamics. A wide range of current sheet widths can be observed. The system presented 
here should be consider as an example of the effect of the boundaries on MHD turbulence. 
To our knowledge, there is a lack of such studies in the literature and further investigation 
along these lines should be extended to other MHD systems (3D MHD, compressible) as 
well. 
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FIGURES 



FIG. 1. Cartoon of the reduced MHD setup. 

FIG. 2. Top: magnetic and kinetic energy (dashed) vs time from a run with S = 2000, 
tA/tp = 0.5. Bottom: Dissipation rate and forcing (dashed) for the same run. 

FIG. 3. Perpendicular energy spectra from S = 2000 runs and different values of the ratio 
t A /t p = 0.1, 0.5, 1. Resolution is 1024 x 1024 x 64, 2048 x 2048 x 64, 2048 x 2048 x 128 respectively. 
On the left panels, magnetic and kinetic energy (dashed) perpendicular spectra. The dotted lines 
show exact power laws given for reference. On the right, total "compensated" perpendicular spectra 
and corresponding power index value a. 

FIG. 4. Spectra and compensated spectra for a run with S = 800, tA/tp = 2. 

FIG. 5. Perpendicular energy spectra for the Elsasser fields z ± from S = 2000 runs and different 
values of the ratio tA/t p = 0.1, 0.5, 1. On the left panels, and (dashed) spectra. On the 
right, compensated spectra for the product and corresponding power index value a. 

FIG. 6. The fc-dependent kinetic Reynolds number = u^jikv) = ^J(E%/k)/i> for a run with 
S = 2000, t A /t p = 1. 

FIG. 7. Resolution study of the spectra for runs with S = 2000, tA/t p = 0.1 

FIG. 8. Cross section of parallel current density for cases with tA/t p = 0.1 and tA/t p = 0.5. 
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